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Executive Summary 


Finite Element Analysis of Reverberation Chambers 

The primary motivating factor behind the initiation of this work was to provide a 
deterministic means of establishing the validity of the statistical methods that are recommended 
for the determination of fields that interact in an avionics system. The application of finite 
element analysis to reverberation chambers is the initial step required to establish a reasonable 
course of inquiry in this particularly data-intensive study. The use of computational 
electromagnetics provides a high degree of control of the “experimental” parameters that can be 
utilized in a simulation of reverberating structures. As the work evolved there were four primary 
focus areas they are 

1 . The eigenvalue problem for the source free problem 

2. The development of a complex efficient eigensolver 

3 . The application of a source for the TE and TM fields for statistical characterization 

4. The examination of shielding effectiveness in a reverberating environment 

One early purpose of this work was to establish the utility of finite element techniques in 
the development of an extended low frequency statistical model for reverberation phenomena. By 
employing finite element techniques, structures of arbitrary complexity can be analyzed due to the 
use of triangular shape functions in the spatial discretization. The effects of both frequency stirring 
and mechanical stirring are presented. It is suggested that for the low frequency operation the 
typical tuner size is inadequate to provide a sufficiently random field and that frequency stirring 
should be used. The results of the finite element analysis of the reverberation chamber illustrate 
the potential utility of a 2D representation for enhancing the basic statistical characteristics of the 
chamber when operating in a low frequency regime. The basic field statistics are verified for 
frequency stirring over a wide range of frequencies. Mechamcal stirring is shown to provide an 
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effective frequency deviation. It can be suggested that a large frequency deviation as compared to 
the operating frequency is required to obtain the desired field statistics. An advantage of 
employing a numerical scheme such as finite elements to analyze reverberation chambers is that 
the field values at virtually any point in the structure can be easily obtained without perturbing the 
field as can occur when using an antenna to monitor the fields. 

Essential to the finite element solution of these problems is an effective numerical 
procedure for solving large-scale, sparse systems of linear equations and generalized eigen- 
equations. These solution phases typically represent the most costly steps of the analysis in terms 
of computational resources. Subspace and Lanczos iterations have been developed, well 
documented, and widely accepted as efficient methods for obtaining p-lowest eigen-pair solutions 
of large-scale, practical engineering problems. The focus of this work is, however, to re-examine 
these popular eigen- solution algorithms, with the viewpoints to incorporate recent developments 
in vectorized sparse technologies in conjunctions with Subspace and Lanczos iterative algorithms 
for computational electromagnetic enhancements. Basic generalized eigen-solution algorithms are 
reviewed. Major computational tasks in subspace iteration and Lanczos algorithms have been 
identified. Efficient sparse technologies have been developed, and fully utilized, in conjunction 
with the basic Subspace iterations and Lanczos algorithms for efficient solutions of the 
generalized eigen-equations. Numerical results from practical finite element models have clearly 
indicated that the proposed “sparse” Subspace iterations and Lanczos algorithms have offered 
substantial computational advantages over the traditional “skyline”, or “variable bandwidth 
strategies. Both lumped and consistent mass formulation have been incorporated into our eigen- 
solution package. 

A two-dimensional finite element model for both the transverse electric (TE) and 
transverse magnetic (TM) solutions inside a simulation reverberation chamber is presented. Tuner 
effects on the modal structure and the resulting statistics of the field distribution have been 
explored. It is demonstrated that the TE simulation reverberation chamber provides excellent 
statistical and reverberation characteristics. The measures of reverberation characteristics 
includes normalized standard deviation, the max-to-average ratio, the stirrer ratio, and the field 



uniformity. In addition the notion of “source” stirring is introduced wherein the fields that were 
coupled into the non-complex geometry were shown to be statistically similar to those in the 
complex environment. The idea of source stirring may shed light on the question of whether it is 
absolutely necessary to stir in all regions for which a reverberation type test is being performed. 
The shielding effectiveness in a reverberating environment was been examined. Coupling 
shielding effectiveness for TE and TM results was presented. The reverberating environment 
provides additional insight into the shielding properties in a statistical sense. 

Shielding of avionics equipment from both interior and exterior electromagnetic threats is 
an important part of the system reliability budget. It is typically assumed that the shielding 
effectiveness of a structure has a level that is a function of fi'equency for various angles of 
incidence. Measurements in an anechoic chamber and analysis using computational methods 
simulate a plane wave environment to quantify shielding effectiveness. An investigation has been 
performed suggesting that an under certain conditions the aircraft fuselage behaves as a 
reverberation chamber when stirred internally. Another investigation revealed that the fields m an 
aircraft flown in the vicinity of a transmitting antenna have levels with characteristics that are 
statistically similar to a reverberation environment. A conclusion that may be drawn from this is 
that an aircraft fuselage is a structure in which the interior fields are not plane waves. This 
significant connection raises an interesting question; What is the shielding effectiveness of an 
aperture in a reverberating environment? Recently published results for the simulation of the 
shielding effectiveness of a two-dimensional TM structure in a plane wave environment. Based 
on the published work, the shielding effectiveness of the TM structure via aperture coupling for 
several frequencies in a simulated reverberation environment are examined. Measurements were 
also performed in NASA’s reverberation chamber and are presented to provide a connection to 
the simulated statistical results. A two-dimensional finite element model for transverse magnetic 
(TM) solution inside a reverberation chamber is presented. Tuner effects on the modal structure 
and the shielding effectiveness in a reverberating environment are exanuned. The shielding 
effectiveness of the TM structure for aperture coupling for several frequencies in a simulated 
reverberation environment are presented and compared to published and measured results. These 
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results point to the development of a statistical model of shielding effectiveness. This model will 
form a critical component to the overall mechanism of upset prediction in digital systems. 

An Example... 

An example of an application of the notion of a probabilistic model for shielding 
effectiveness a problem of determining the probability of exceeding a particular voltage for a 
given input voltage. Consider the problem of determining the probability of having an electric 
field intensity of 1 V/m anywhere inside a cavity exposed at 1000 V/m. The shielding effectiveness 
has been characterized as a Normal random variable with mean of 35.3 dB and standard deviation 
of 10.5 dB at a frequency of 400 MHz. The probability of exceeding IV/m for a 1000 V/m 
exposure corresponds to calculating the probability of having a shielding effectiveness less than 60 

dB which is F (60) = 1 -of = 0.9906 . This calculation predicts a 99.06% probability 

SE\ ) JQ5 j 

of exceeding IV/m for an exposure of 1000 V/m. If the threat level is reduced to 100 V/m the 
probability of exceeding 1 V/m reduces to 67.36 %. 
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Chapter 1 

A Two Dimensional Finite Element Analysis of 
Reverberation Chambers 

1.1 INTRODUCTION 

A reverberation chamber is an enclosure consisting of metal walls with a metallic paddle 
wheel (denoted a “stirrer” or “tuner”) essentially forming a high quality factor (0 cavity with 
continuously variable boundary conditions. The fields inside the chamber for a given tuner 
position are completely deterministic. However, given the nature of the variable boundary 
condition, the ability of a given source to couple energy into certain modes, and the passband 
characteristic due the chamber 0 the fields are typically characterized by statistical means. 
Specifically the probability density functions for the real and imaginary parts of the electric and 
magnetic fields are normally distributed [1]. In a general three dimensional cavity the field 
magnitude is j^with six degrees of freedom - the real and imaginary parts of the three field 
components. A model for the analysis of a two-dimensional cavity consists of an examination of 
a waveguide structure evaluated at the cutoff frequencies for the modes. Note that at the cutoff 
frequency the mode exists but does not propagate in a manner similar to the cavity modes. For 
this study at cutoff the field magnitudes are Rayleigh distributed and the squared fields (related 
to the power) are exponentially distributed [1] since at most two degrees of freedom exist. The 
conformity of the actual fields to the specific statistical description is dependent upon the number 
of modes that can simultaneously exist at a given frequency. The number of modes that can exist 
at a given frequency is a function of the cavity Q and the operating frequency. Hill [3] gives the 
two-dimensional specific mode density for finite Q as 

^ ^ 2 ^^ ( 1 . 1 ) 

* Qc^ 

with a and Z) the dimensions of the structure,/ is the operating frequency, c is the velocity of 
light. 



Although the primary methods of characterizing reverberation chambers have been 
experimental there have been at least two papers which address the numerical analysis of two- 
dimensional reverberation chambers. The transmission line matrix (TLM) method was 
employed by Wu [2] to examine the characteristics of the fields in the presence of a large stirrer 
in a reverberation chamber. Wu used a two-dimensional cavity with a one-dimensional stirrer 
model to examine the shifting of the eigenfreq uencie s and the effects on the mode amplitudes for 
a variety of different stirrer sizes. Wu concludes that the frequency shift for a given stirrer 
essentially determines the random nature of the fields. Wu also makes an important observation 
regarding the consequences of the shifting eigenfrequencies in noting that the effect of a large 
stirrer is essentially a form of random modulation. This modulation is considered to have both 
amplitude and frequency modulation effects. Wu concludes that an ineffective (electrically 
small) stirrer is unable to provide adequate frequency modulation effects. Hill [3] examined the 
concept of electronic mode stirring initially suggested by Loughry [4]. Hill employs a Green s 
function approach to the analysis of an empty two-dimensional rectangular cavity with an 
electric line source. Electronic mode stirring essentially implies that instead of stirring via 
mechanical means to obtain a statistical field distribution, the source can be swept over a narrow 
bandwidth to effect the shift in eigenfrequencies as suggested by Wu [2]. 

The purpose of this paper is to establish the utility of finite element techniques in the 
development of an extended low frequency statistical model for reverberation phenomena. By 
employing finite element techniques, structures of arbitrary complexity can be analyzed due to 
the use of triangular shape functions in the spatial discretization. The TLM method uses a 
rectangular grid structure which introduces an error in the discretization of structures which do 
not lie along a Cartesian coordinate axis. The finite element method for the source free two- 
dimensional cavity model yields the fields at all discretized points within the geometry as a result 
of solving an eigenvalue problem. By contrast the TLM method requires a Fourier transform of 
the data at all points desired in the field representation. Hill’s analysis is applicable to an empty 
structure so that the relationship between mechanical stirring and frequency stirring cannot be 
readily explored. Both Wu and Hill consider a two-dimensional structure with a source at a 
sufficiently high frequency to assure random fields and reasonable field homogeneity consistent 
with an overmoded structure. The field statistics for reverberation chamber are well known for 
the case when the chamber is clearly overmoded (such as when exciting the chamber at with an 
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operating frequency that is six times higher than the first resonant frequency of the chamber), but 
are less defined for undermoded low frequency chamber operation. 

In the following section the basic finite element scheme will be described. The approach 
uses edge elements for the transverse field to eliminate spurious solutions. In section 3 the 
results for a rectangular structure will be examined. The effects of both frequency stirring and 
mechanical stirring will be presented. It can will be suggested that for the low frequency 
operation the typical tuner size is inadequate to provide a sufficiently random field and that 
frequency stirring should be used. The strengths and limitations of the approach will be 
presented in the conclusion. 

1.2 FINITE ELEMENTS FOR THE REVERBERATION CHAMBER [5-7] 

A model for the analysis of a two-dimensional cavity in this paper will consist of a 
waveguide structure evaluated at the cutoff frequencies for the modes. Note that at the cutoff 
frequency the mode exists but does not propagate. In the two-dimensional model to be used the 
cutoff frequencies can me made directly analogous to the resonant frequency of the three 
dimensional structure. One significant distinction is that the cavity can be considered to possess 
a bandpass response centered at the operating frequency provided that a cavity mode is supported 
at that frequency. The width of the response in inversely proportional to the quality factor (0 of 
the structure. In contrast the two-dimensiona l w aveguide structure can be considered to have a 
high pass response which supports all modes with cutoff frequencies less than the operating 
frequency. By using a lossless structure it is possible to examine the infinite Q case where only 
one mode is considered to exist. It is suggested that this model adequately reflects the physical 
situation for the low frequency reverberation chamber. 

For the development of the finite element method Maxwell's curl equations can be 
considered in a source free, homogeneous, isotropic, time harmonic form (all sources are 

assumed steady-state with time dependence) as given by 
VxE = -jco^H , 
and WxH = jco£E. 

Although analytical results are available for a large class of problems for static, quasi-static, and 
dynamic conditions, numerical methods typically must be applied whenever the geometry does 
not coincide with a separable coordinate system. 


3 



A full-field approach is required when the medium is characterized by higher order 
modes, hybrid modes such as for a general inhomogeneous media, or a waveguide or cavity with 
imperfectly conducting walls. Other problems can be expressed in a scalar form when there is a 
distinct separation between TE and TM modes. In general this approximation cannot be applied 
and a full vector formulation is required. The solution for the fields inside the reverberation 
chamber will be required to support hybrid modes in the presence of any metallic discontinuity 
or an inhomogeneity due to lossy objects. For the finite element approach the electric field in 
terms of all three vector components is the desired field quantity. A two dimensional 
representation for the fields may be expanded in terms of the transverse field and the longitudinal 
field where 

E (x, y, z) = [e, (x, y ) + ze^ (x, y)] . (1.4) 

Expressing the del operator as V = V, + and setting y = 0 for cutoff, it is possible to write 


two separate equations - one for the transverse part and another for the z component in terms of 
the cutoff frequency. 


V.x 


— V,x?, 


}^r 


-kle^e,=0, 


V.e, 


— k^ee = 0 , 

c r z ^ ? 


leading to the following weighted residual form: 

JJ— (V, xe,).(v, x7;)t/5 = .7; c/5 

Mr 


(1.5) 

( 1 . 6 ) 

(1.7) 


and 


Jf— (>-8) 

Where = 2;r and is the known as the cutoff wavenumber. Setting y = 0 in eq. (1.4) has 

the effect of rendering the general three-dimensional problem into a two-dimensional problem by 
suppressing the z-behavior. These modes do not propagate and are analogous to the excitation 
frequencies of a two-dimensional cavity. 

The transverse components have been discretized using edge elements and the 
longitudinal component uses node-based Lagrange functions. After integrating and summing 
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over all elements in the domain the following eigenvalue problem of the form Ax = ABx for the 
cutoff wavenumber is obtained 


5, 0 

0 5 , 



0 



(1.9) 


The solution to the matrix eigenvalue problem yields an eigenvalue-eigenvector pair 
where the eigenvalue represents the cutoff wavenumber, , and thereby the cutoff frequency. 


/^, for the corresponding eigenvector. The eigenvector represents the field solution that is 

supported at frequencies at and above the cutoff frequency. The reader should be aware that by 
obtaining the solution to an eigenvalue problem it is implied that for a given excitation frequency 
it is assumed that the corresponding mode will exist . The basic nature of this source free 
problem is that a mode exists by virtue of applying a frequency above the cutoff frequency. In a 
structure with a source only those modes that are supported by the source will exist. 
Additionally, since we are at cutoff the implication is that each mode exists independently. This 
model corresponds to a cavity with an infinite Q. The following section will present the results 
of the finite element analysis to the lossless two-dimensional reverberation chamber model. 


1.3 RESULTS 

This paper will consider the basic two-dimensional geometry as depicted in Figure 1.1. 
An anti-symmetric tuner may be rotated to provide mechanical mode stirring, however, the 
primary results of this work examine the use of a frequency stirring approach. A FORTRAN 



Figure 1.1. Basic 2D Reverberation Chamber. 
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code applying finite element techniques has been written that provides the electric field at any 
point within the chamber. 

The challenge of applying finite element techniques to reverberation chambers is that 
enormous computational resources are required. The primary limitation is in the available 
random access memory (RAM) required for in-core processing. Virtual storage, or out-of-core 
processing, of the data can also be used, but greatly slows the solution time. Taking advantage 
of sparsity in the resulting system of equations helps significantly. In spite of the advantages 
gained by using the sparse eigenvalue solver optimized for the CONVEX the system of 
equations was limited to the value used in this paper at around 1 0,000x 1 0,000. (The CONVEX 
220 at NASA-Langley has 1 GB of random access memory.) By taking advantage of sparsity 
the dimensions of the matrices A and B in eq. (1.9) are reduced to approximately 10,000x7. 
However, the computation of the eigenvectors results in a full solution matrix whose dimensions 
are dependent upon the number of eigenvectors desired. 

For the reverberation chamber simulation a two-dimensional model of National 
Aeronautics and Space Administration (NASA) Langley Research Center's B chamber with 
dimension 3.96 m x 7.10 m was generated and is depicted in Figure 1.2. The resulting system of 
equations was on the order of 10,000x10,000 and was solved using a CONVEX 220 using a 
sparse eigenvalue solver in FORTRAN. The data was then processed in MATLAB. The 
solution of eq. (1.9) has on the order of 10,000 eigenvalues and eigenvectors as possible values 
of cutoff frequency and corresponding field configurations. Although there are approximately 
10,000 solutions available, only those solutions with cutoff frequencies that permit a 
discretization of 10 samples per wavelength were retained. The rules governing sampling 
density are not strict and can vary from problem to problem [8] md are related to the Nyquist 
rate for sampled data. The sampling density has small effect on the accuracy of the cutoff 
frequency, but has a significant effect on the accuracy of the corresponding field distribution. 
The level of discretization allowed accurate solutions up to roughly 300 MHz, and resulted in a 
possible 171 valid cutoff frequencies. A typical result for the field behavior is shown in Figure 
1.3 for the x-component of the electric field, at a frequency of 108.1 143 MHz. Note that the 
highest field levels are in the immediate vicinity of the stirrer as would be expected since the 
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fields become infinite near the tuner. Figure 1.4 depicts a higher frequency, 301.069 MHz, for 
the infinite Q structure and depicts a significantly more complex field structure. 



Figure 1.2. Discretization for the two-dimensional geometry. 
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Typical experimental work using a reverberation chamber involves the rotation of the 
tuner and a measurement of the received power to determine such characteristics as the stirring 
ratio, or a measurement of the voltage standing wave ratio (VSWR) based on reflected power at 
the transmit port. The rotation of the tuner provides a variable boundary condition that results in 
a deterministic variation of the field by applying finite elements. One issue of concern for the 
eigenvalue problem of eq. (1.9) involving the cutoff frequency is that although all modes are 
found corresponding to the particular geometrical configuration, the weighting of the modes is 
not possible. In an actual reverberation chamber the source does not couple energy equally 
between the possible modes. Another concern is that the modal fields represented by the 
eigenvectors exist at discrete frequencies and the chamber Q does not directly enter into the 
formulation. The implication is that the computed modal fields are for an infinite Q chamber, 
and that multi-moding strictly speaking does not occur since only one mode can exist at a time. 
We can consider Q = fojB for a passband structure with /o corresponding to the cutoff 
frequency fc- An infinite Q would imply that the bandwidth B would be zero, thus allowing only 
one mode to exist. 

1.3.1 Frequency Stirring 

Motivated by the notion that the action of stirring excites the modes weighted by a factor 
governed by the quality factor of the chamber, and noting that frequency stirring may be used in 
place of mechanical stirring [3] an interesting experiment was performed. For a fixed tuner 
position computing all the eigenvalue-eigenvector pairs, within the accuracy constraint based on 
proper discretization, corresponds to computing all the possible modes that could he excited in 
the chamber. The solutions for cutoff result in a field that is either transverse electric (TE) or 
transverse magnetic (TM) as the frequency is varied. Consequently, when the solution is TE 
only Ex and Ey exist and when the solution is TM only E^ exists. Figure 1 .5 depicts the amplitude 
oi Ex at one point in the chamber as the frequency is varied. The field variations increase with 
frequency. This is related to the increased mode density at higher frequencies. A histogram was 
constructed using MATLAB in Figure 1 .6 based on the data of Figure 1 .5. 
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- 0.5 > 0.4 - 0.3 - 0.2 - 0.1 0 0.1 0.2 0.3 0.4 0.5 

E (TE) amplitude, 20 bins, xbar =-0.83694, ybar =-l .61 7 


Figure 1.6. The histogram of the data of Figure 1.5. 



Any single component should have a Gaussian, or normal distribution [1] with the probability 
density function, or pdf, given by 




( 1 . 10 ) 


with mean // and standard deviation a. The results obtained from the histogram are 
approximately Gaussian considering the limited number of data points used (101 total points for 
TE). The histogram provides one means of checking basic trends in the data as compared to a 
known probability density function, but is limited in its ability to adequately characterize the 
distribution. A useful visualization tool is the use of a probability plot. The normal probability 
[9,10] plot essentially renders the cumulative distribution of a normal distribution to that of a 
straight line as depicted in Figure 1.7 (a) for a set of normally distributed data generated from a 
random number generator. The data are plotted with “+” symbol. A solid line connects the first 
and third quartiles of the data corresponding to a robust linear fit of the sample statistics. The 
linear fit is robust with regard to extremes in the data. The dashed line on the plot represents an 
extrapolation of the data to the ends of the sainple. Data that comes from a normal distribution 
will appear linear on the plot. If the sample data comes from another distribution the data appear 
to add curvature in the plot. Consider the curvature in Figure 1.7 (b) for a set of exponentially 
distributed random data. The curvature in Figure 1 .7 (b) clearly demonstrates that the data are 
not normally distributed. The x-component of the electric field, fi, is examined by means of a 
probability plot in Figure 1.8. The linearity observed is sufficient to assert that the is normal 
for the purposes of applying procedures that assume that Ex is normal, such as combining with 
other field components to assess the statistics of the squared total field as will be done later. A 
similar test is performed on the other TE field component Ey in Figure 1 .9. For the TM fields Ez 
is obtained and is validated in Figure 1.10. 
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Figure 1.7. Normal probability plot from a random number 
generator for a data count comparable to the data obtainedfor the 
fields for (a) normal data, and (b) exponentially distributed data. 
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Normal Probrf>ilil> Plot. E^(TE) 




Dm. =-0 0081204 0^=0.022565, o =0.15022 

Figure 1.8. Normal probability plot for Ex for a single location as the 
frequency is varied. 
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Normal Probability Plot, E^, (TE) 



Data. h=- 0.018494, =0,020379. o=0.14276 


Figure 1.9. Normal probability' plot for E, for a single location 
as the frequency is varied. 
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Normal Probability Plot, (TM) 



Data, M =0.02408, a' =0.046052, a =0.2146 


Figure 1.10. Normal probability plot for for a single location as 
the frequency is varied. 


Parameter estimation can be used as a further validation of the underlying statistics by 
determining the maximum likelihood estimators (MLEs) for data assumed to be taken from a 
particular distribution. For a normal distribution the MLEs are the mean and standard deviation. 
The 95 % confidence intervals for the mean Md variance are also of interest. The MLEs of a 
sample set should lie within the confidence interval to assert reasonable agreement of a sample 
distribution with a known probability density function. The results of the parameter estimation 
are detailed in Table 1 showing the MLEs and confidence intervals for the data of £v, Ey, and E^. 
The parameter estimation on the field data generated as a result of frequency stirring 
demonstrates significant agreement within the 95 % confidence interval. The MATLAB 
Statistical Toolbox was used for the parameter estimation.’ 
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Table 1. Maximum likelihood estimators (MLEs) and 95% confidence intervals. 


Normal 

f^St 

^ki 


^est 

E. 

-0.00812 

(-0.37780,0.02150) 

0.15022 

(0.13197,0.17436) 

Ey 

-0.01849 

(-0.04668,0.00969) 

0.14276 

(0.12542,0.16570) 

E. 

0.02408 1 

(-0.02709,0.07529) 

0.21460 

(0.18400,0.25750) 

Exponential 

^est 

dci 

best 

bci 

El+El 

32.381 

(25.763,38.998) 

1.1870 

(-14.523,16.761) 


Based on the results of the normal probability plot it can be asserted that both and Ey 
are normal and can combined as the sum of squared normal data which will result in an 
exponential distribution. For the TE fields this squared field (proportional to power) data is 

c = £.;+£; ('■'» 

The squared field data is depicted in Figure 1.11 with the corresponding histogram in Figure 
1.12, and resembles the expected exponential distribution [1]. 


(TE). Element Number =5021 



Figure 1.11. The squared field (power) for a single point in the reverberation 

chamber. 
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Histogram for (TE), Element Number -502 1 

25 1 ^ ^ ^ ^ ^ ' ' ' ' 



E^(TE) amplitude, 20 bins, xbar =-0.83694, ybar — 1 .61 7 
Figure 1.12. The histogram for Figure 1.11. 


The probability density function for an exponential distribution is 
y = f{x\^i) = ^e ^ 

with parameter //. In order to assess the agreement of the total squared field data a Weibull 
probability plot can used in exactly the same manner as the normal probability plot. The 
exponential distribution is a special case of a Weibull distribution given by the pdf 

= ,x>0 (113) 

■^1 0 , x <0 


E 

L 
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The exponential distribution corresponds to 6 = 1 in eq. (1.13). The squared field data are 
analyzed by means of the Weibull probability plot in Figure 1.13 and demonstrates good 
agreement. The MLEs for a Weibull distribution are the values of a„, = 32.3807 and bes, = 1 • 1 87 
in eq. (1.13). Of particular significance is the estimator for b which is relatively close to that of 
the exponential distribution. The 95 % confidence intervals for the estimates are provided in 
Table 1. 


Weibull Probability Plot, (TE) 



Figure 1.13. Weibull probability plot for the squared fields for a single location as 

the frequency is varied. 


The effects of mode density can be seen in Figure 1.14 for the change in frequency. A/, 
versus frequency. Note that the mode density increases as the frequency is increased. One 
implication of increasing mode density is that additional modes are supported for a fixed non- 
zero bandwidth giving rise to a multimode environment. For an infinite Q structure the 
bandwidth is essentially zero so that this implication cannot be readily explored under the current 
source free environment. Another implication increased mode density is that the field can be 
significantly altered by exciting alternative modes via changing the cutoff/resonant frequencies 
of supported modes by changing the physical structure, thus altering the boundary conditions. 
The changing of boundary conditions is typically accomplished by mechanical stirring and will 
be discussed in the next subsection. 
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Figure 1.14. A plot of Af vs. f illustrating the mode density increase as frequency 

increases. 


1.3.2 Mechanical Stirring 

The action of mechanical stirring is an important aspect of reverberation chamber 
analysis. Applying finite elements to a rotating stirrer requires a significant amount of data 
processing since the field will need to be computed for each tuner position. For the purposes of 
this investigation the tuner was rotated in 1. 6° increments for a total of 225 steps. The increment 
was chosen to match the step size used by the National Institutes of Science and Technology 
(NIST) when NASA Langley’s chambers were characterized in early 1997. The finite element 
discretization was reduced to a system size of 3,000 unknowns resulting in an upper frequency 
limit of 150 MHz. The eigenvalue problem of eq. (1.9) was solved for each position of the tuner 
which resulted in a total of 43 eigenvalues between the frequencies of 20 MHz and 150 MHz. 
Figure 1.15 (a) depicts the total field at approximately 60 MHz for three independent points 
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located far from the tuner and greater than half a wavelength from any wall. The total field is 
periodic and relatively smooth. Figure 1.15 (b) shows similar results for a frequency of 
approximately 47 MHz. 

Work was undertaken to explore various effects of the tuner on the frequency difference 
between each mode. In this work the assumption of cutoff is advantageous since it is possible to 
look at each mode in isolation. The first study is the change in frequency versus mode number 
for a fixed tuner position. The results of Figure 1.16 (a) depict a plot of the frequency difference 
between the modes for a 1.6° shift in the tuner at a final tuner position of 32 degrees 
(representing the frequency difference as the tuner is moved from 30.4 to 32 degrees - an 
arbitrary selection) in part (a). Figure 1.16 (b) depicts the absolute peak frequency deviation for 
all tuner positions. This study may provide insight into the degree to which a particular mode is 
affected by the presence of the tuner. One question is; Is there a mode that is more greatly 
affected by the tuner than another? Another question is; Which mode exhibits the largest 
maximum frequency deviation? It can be suggested from the data of Figure 1.16 (b) that modes 
19 and 20 are most greatly affected by the tuner. The results are suggestive of a means of 
addressing a long standing question [11] regarding the amount of frequency shift of a given 
mode that can be obtained by stirrer action to facilitate a comparison to frequency stirring. 

The notion of frequency modulation effects [2] suggest a second study to analyze the 
frequency difference for a particular mode versus tuner position. This second study resulted in 
the data as presented in Figure 1.17. It can be seen that the modes are significantly affected by 
the tuner in a complex manner that similar to that of frequency modulation. Figure 1.17 (a) 
presents the frequency deviation of a single mode as a function of tuner position. Figure 1.17 (b) 
depicts the absolute peak frequency deviation for all modes versus tuner position. From Figure 
l.s 16 and 17 it can be suggested that the tuner exhibits a peak instantaneous frequency deviation 
of about 1 .4 MHz. The random deviation is desired to excite as many modes as possible at a 
given frequency. 

As a third study, animations of the field have been made in which each frame 
corresponds to the solution for a particular tuner position. These animations provide a means of 
checking mode continuity for the fields obtained from the matrix of eigenvectors. The modes 
can shift as the tuner is rotated and it is useful to examine those modes for which a radical 
departure from a particular mode shape occurs. For some modes the fields rapidly change from 
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TE to TM as the tuner is rotated. This mode shifting occurs most significantly for modes 19 and 
20 as can be inferred from Figure 1.16. This maximum frequency deviation modes 19 and 20 is 
nearly 1 .4 MHz. This large deviation is indicative of mode sensitivity due to the tuner position. 
It can be suggested that this is a desired condition for proper reverberation chamber operation in 
a multimode environment. As the cutoff frequency shifts around a particular operating 
frequency, modes are excited or attenuated (they become evanescent in the two-dimensional 
structure) and the field variation provides the desired overall field statistics and homogeneity. 
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Figure 1.15. A plot of the total field for 225 tuner steps at two different points in the 


chamber, (a) The total field at approximately 60 MHz for three independent points 


located far from the tuner and greater than half a wavelength from any wall, (b) 


Similar results at a frequency of approximately 47 MHz. 
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Figure 1.16. Frequency deviation versus the modes, (a) One tuner position, (b) All 

tuner positions. 
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Figure 1.17. Frequency deviation versus tuner position, (a) Single 


mode (b) All modes. 
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1.4 CONCLUSIONS AND FUTURE WORK 

This paper has presented a two-dimensional analysis of reverberation chambers using 
finite element techniques. The results of the finite element analysis of the reverberation chamber 
illustrate the potential utility of a 2D representation for enhancing the basic statistical 
characteristics of the chamber when operating in a low frequency regime. The basic field 
statistics are verified for frequency stirring over a wide range of frequencies. Mechanical stirring 
was shown to provide an effective frequency deviation. It can be suggested that a large 
frequency deviation as compared to the operating frequency is required to obtain the desired field 
statistics. The field statistics of mechanical stirring was not specifically explored in this paper 
since the modes for the cutoff exist in isolation and do not yield desirable characteristics. An 
advantage of employing a numerical scheme such as finite elements to analyze reverberation 
chambers is that the field values at virtually any point in the structure can be easily obtained 
without perturbing the field as can occur when using an antenna to monitor the fields. 

Work is continuing in the development of a high speed sparse complex eigensystem 
solver that will allow the solution of much larger systems (with on the order of 60,000 
unknowns) to be analyzed. Other current work includes the coupling of a source to the fields in 
the 2D geometry by solving the eigenvalue problem for the propagation constant and including a 
source. The resulting solution would contain the proper mode weighting for the particular 
configuration depending upon the stirrer position. The author is currently examining the two- 
dimensional line source problem analogous to Hill’s [3] work that addresses low frequency 
effects on field homogeneity. Future work will also examine the three dimensional problem at 
low frequencies that include the source. 
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Chapter 2 

Subspace and Lanczos Sparse Eigen-Solvers for Finite 
Element Structural and Electromagnetic Applications 


2.1. INTRODUCTION 

The finite element method has been used successfully for the solution of many practical 
engineering problems in various disciplines, such as structural analysis, fluid mechanics, structural 
optimization, heat transfer, electromagnetic etc. [1-7], Essential to the finite element solution of 
these problems is an effective numerical procedure for solving large-scale, sparse systems of linear 
equations and generalized eigen-equations. These solution phases typically represent the most costly 
steps of the analysis in terms of computational resources. 

Subspace and Lanczos iterations have been developed, well documented, and widely 
accepted as efficient methods for obtaining p-lowest eigen-pair solutions of large-scale, practical 
engineering problems [1,2]. The focus of this paper is, however, to re-examine these 2 popular 
eigen-solution algorithms, with the viewpoints to incorporate recent developments in vectorized 
sparse technologies in conjunctions with Subspace and Lanczos iterative algorithms for 
computational electromagnetic enhancements. 

A reverberation chamber (see Figure 1) is an enclosure consisting of metal walls with a 
metallic paddle wheel (denoted a “stirrer” or “tuner”) essentially forming a high quality factor (Q) 
cavity with continuously variable boundary conditions. The fields inside the chamber for a given 
tuner position are completely deterministic. However, given the nature of the variable boundary 
condition, the ability of a given source to couple energy into certain modes, and the passband 
characteristic due the chamber Q, the fields are typically characterized by statistical means. 
Specifically the probability density functions for the in-phase and quadrature component of the 
electric and magnetic fields are normally distributed [8]. The field magnitudes are Rayleigh 
distributed and the power is exponentially distributed [8]. The degree to which the actual fields 
conform to the specific statistical description is dependent upon the number of modes that can 
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simultaneously exist at a given frequency is a function of the cavity Q and the operating frequency. 

By utilizing a numerical analysis scheme such as finite elements the effects of introducing 
loss maybe accurately simulated, however, for this paper a two-dimensional lossless simulation will 
be performed. 



Figure 1. Basic 2D Reverberation Chamber. 

A two-dimensional model may yield insight into the nature of the loss experienced on real airframes 
and also provide an error measure on data in the ideal mode-stirred chamber. 

2.2 FINITE ELEMENTS FOR THE REVERBERATION 
CHAMBER [7] 

Consider Maxwell’s curl equations in a source free, homogeneous, isotropic, 
form given by 

VxE = -jco^y 

V xH = j6)£E, 


time harmonic 

( 2 . 1 ) 

( 2 . 2 ) 
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Although analytical results are available for a large class of problems for static, quasi-static, and 
dynamic conditions, numerical methods typically must be applied whenever the geometry does not 
coincide with a separable coordinate system. 

A fiill-field approach is required when the medium is characterized by higher order modes, 
hybrid modes such as for a general inhomogenwus media, or a waveguide or cavity with imperfectly 
conducting walls. Other problems can be expressed in a scalar form when there is a distinct 
separation between TE and TM modes. In general this approximation cannot be applied and a full 
vector formulation is required. The solution for the fields inside the reverberation chamber will be 
required to support hybrid modes in the presence of any metallic discontinuity or an in homogeneity 
due lossy objects. For the finite element approach the electric field in terms of all three vector 
components is the desired field quantity. A two dimensional representation for the fields may be 
expanded in terms of the transverse field and the longitudinal field where 


E[x,y,z)= e,{x,y) + ze^(x,y) 




( 2 . 3 ) 


Expressing the del operator as V = and setting - 0 for cutoff, it is possible 

to write two separate equations, one for the transverse part and another for the z component in terms 
of the cutoff frequency, 

V,x(^V,xe,)-*,V, = 0 (2'') 

leading to the following weighted residual form: 

niil7(v,xe,)»(v,x7;)z/5= k^e^jje,»T,ds (2.6) 

and 

( 2 . 7 ) 
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The transverse components will be discretized using edge elements and the longitudinal 
component used node-based Lagrange functions. After integrating and summing over all elements 
in the domain, the following eigenvalue problem of the form Ax = >, Bx for the cutoff frequency 
is obtained 


1 

o 

1 


II 

'T, 

1 
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[0 


€ 

0 


— 1 
1 


( 2 . 8 ) 


2.3. BASIC SUBSPACE ITERATION ALGORITHM [1,3] 


The generalized eigen-equations (shown in Eq. 2.8), in matrix notation, can be expressed as 


lK][<p]=lM][4X\ 


(2.9) 


In Eq. (9), matrices IX] and |>l]r^)resait die stmcturalstiffiiess and nrass,te^)ective Matrices [A] and [^] 


represent the eigenvalues and eigenvectors, respectively. The dimension (or degree-of-freedom) of 
matrices in Eq. (9) is N. For many practical engineering applications, [K] is symmetrical and 
positive definite. Subspace iteration algorithm can be used effectively to obtain the lowest p eigen- 
pair solutions. The algorithm can be conveniently described by the following step-by-step 
procedures as shown in Table 1 . 
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Table 1: Step-by step Basic Subspace Algorithm 


Step 1 ; Select the starting iteration vectors [Y,] where q « N 
Step 2; Factorize the structural stiffness matrix 

[K] = [L][D][LY (2.10) 

In Eq. (1 0), [ L ] is the lower triangular matrix, and [ D ] is the diagonal matrix 

Step 3: For k = 1 ,2, Maxiter, where Maxiter represents the input maximum number 

of iterations, the following tasks need to be done 

Step 4: Solve ] ^xq froiri the following matrix equations 



(2.11) 

Step 5: Compute the reduced stiffness matrix 

(2.12) 

Step 6: Compute the reduced mass matrix 



(2.13) 


(2.14) 

Step 7: Solve the reduced eigen-equations 

(2.15) 

The eigenvalues and the associated eigenvectors [ Q k+i ] need to be arranged 

in the ascending orders (for example ••• ) 

Step 8; Find an improved approximation to the eigenvectors 

(2.16) 


Step 9: Check for convergence. The iterative process will be stopped if either 
convergence is achieved, or the maximum number of iteration ( = Maxiter) is reached 
(or else, return back to step 3). 


30 



2.4. LANCZOS ALGORITHM [1,2] 


Recently, the Lanczos algorithm for the solution of generalized eigenvalue problems has been 
receiving a lot of attention due to its computational efficiency . The original, generalized eigenvalue 
equation can be written as: 

AT (p = (O^Af(p (2.17) 

or 

= ®>(p (2.18) 

where K and M are structural stiffness matrix and mass matrix, respectively, K„= K - oM, a is the 
shift value and =0)^ - a 


Instead of solving Eq. (2.17), or Eq. (2.18) directly, the Lanczos algorithm generates a tri-diagonal 
matrix T„, 


a, P2 

P2 “2 P3 

P3 “3 


f'm 

P„ 


(2.19) 


through the following three-term recurrence: 

r. = 

or in the matrix form: 

[K;'M\Q„ - QJ, = {0.0....r,) = r„e\ 
T z = ^ 

m 


( 2 . 20 ) 

( 2 . 21 ) 

( 2 . 22 ) 


where e^„, = (0,0,... 1 ), Q„, is a Nxm orthogonal matrix with columns qj = 1 ,2,3 ...m, and m is usually 
much smaller than N. By solving the following reduced eigen system, the eigen solution of Eq. 
(2. 1 8) can be obtained as 
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03 . 


e 


9 ^ Qn? 


(2.23) 

(2.24) 


For most structural engineering problems, only a few lowest frequencies and the 
corresponding mode shapes are required, so we have m « N which leads to a significant savings 
in the number of operations. 

A partial restoring orthogonality scheme and a convergence criterion are developed and 
incorporated into the basic Lanczos algorithm, which is described in a step-by-step procedure, shown 
in Table 2. 


Step 1. 

Step 2. 
Step 3. 


Step 4. 


Step 5. 
Step 6. 
Step 7. 
Step 8. 
Step 9. 
Step 10. 


Step 1 1 . 


Step 12. 
Step 13. 


Table 2: Step-by-Step Basic Lanczos Algorithm 

Factorization : = A A A ^ 

Form starting vector: Yo^ 0; Oq = 0 
Compute: M yo 

Compute : 

P. = ^ 

Compute : n, = MG, 

Lanczos iteration 
Forj = 1,2,3, . . . ,do 

^ip = 

\ = e,-P^6^, 

= 0/ M 6, = 

Y.P = 5, - 
= M Y.P 

Reorthogonalization of qj+, 


IF necessary solve Eq( 14) : T^^ = 0^ 
Converged? ( If "No", then return to step 5) 
Eigenvector transformation: <{) = 0^^ 
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2.5. COMPUTATIONAL ENHANCEMENTS FOR SUBSPACE AND LANCZOS 
ALGORITHMS 

It has been pointed out in Sections 2.3-2A that matrix factorization, forward & backward 
equation solution, and matrix-vector (or matrix-matrix) multiplications represent the major 
computational lasts for Subspace iteration, and Lanczos algorithms. Recent developments in 
Sparse technologies [3] will be fully utilized to improve the computational efficiency of both 
Subspace iteration, and Lanczos algorithms. 

2.6. BASIC EQUATION SOLUTION ALGORITHMS 

The Choleski (or U^U) factorization is efficient, however, its application is limited to the 
case where the coefficient stiffness matrix [K] is symmetrical and positive definite. With 
negligible additional computational efforts, the LDL^ algorithm can be used for broader 
applications (where the coefficient matrix can be either positive, or negative definite). In this 
algorithm, the given matrix [K] in Eq. (1) can be factorized as 

[K]=[l][d][lY (2.25) 

where [L] and [D] are lower triangular matrix (with unit values on the diagonal), and diagonal 
matrix, respectively. For a simple 3x3 symmetrical stiffness matrix, Eq. (2.25) can be explicitly 
expressed as 



K ^2 



1 

0 

0 


D, 

0 

0 


1 


^31 ' 

^2. 

K22 

^23 

= 

^2. 

1 

0 


0 

P2 

0 


0 

1 

^32 

.^ 3 . 

Kn 

^ 33 . 


.^ 3 . 

L32 

1 


0 

0 

^ 3 . 


0 

0 

1 


The unknown Lj^ and D| can be easily obtained by expressing the equalities between the upper 
triangular matrix (on the left-hand-side) and its corresponding terms on the right-hand-side of Eq. 
(2.26). 


33 


Since the LDL’^ algorithm will be used later on to develop efficient, vectorized sparse algorithm, 
a pseudo-FORTRAN skeleton code is given in Table 3A (assuming the original given matrix is 
symmetrical and full). 

1 . c . . . Assuming row 1 has been factorized earlier 

2. DO 1 1 I =2, N 

3. DO 22 K =1,1-1 

4. c . . . Compute the multiplier (Note: U represents L^) 

5. XMULT =U(K,I)/U(K,K) 

6. D0 33J=I,N 

7. U (I, J) = U (I, J) - XMULT * U (K, J) 

8. 33 CONTINUE 

9. U (K, I) = XMULT 

10. 22 CONTINUE 

11. 11 CONTINUE 

Table 3A: Skeleton FORTRAN Code for LDL^ 

(assuming the matrix U is completely full) 

2.7. STORAGE SCHEMES FOR THE COEFFICIENT STIFFNESS MATRIX 

Successful implementation of a sparse equation solution algorithm depends rather heavily 
on the reordering method used. While the Reversed Cuthill-Mckee (RCM), or Gipspoole- 
Stockmyer (GS) . . . reordering algorithms can be used efficiently in conjunction with skyline or 
variable bandwidth equation solution algorithms [9-10], these reordering algorithms are not 
suitable for sparse equation solution algorithms. Designing efficient sparse-reordering 
algorithms is a big task itself, and is outside the scope of this paper. For complete treatments on 
this subject, the readers are strongly recommended to popular textbooks and articles in the 
literature [9-10]. In this section, it is assumed that the best available sparse-reordering algorithm. 
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such as Modified Minimum Degree (MMD), or Nested Dissection (ND) [13], has already been 
applied to the original coefficient matrix [K]. To facilitate the discussions in this section, 
assuming the resulted matrix [K] (after using MMD, or ND algorithm) takes the following form 


(2.27) 


In the sparse storage scheme, only non-zero values (and their corresponding locations) of the off- 
diagonal terms need to be stored. The diagonal terms will be stored in a separated array. For the 
data shown in Eq. 27 the number of non-zero off diagonal terms (NCOEF) is 6 (thus NCOEF=6). 


2.8. SPARSE SYMBOLIC FACTORIZATION 

The purpose of symbolic factorization is to find the locations of all nonzero (including 
“fills-in” terms), off-diagonal terms of the factorized matrix [U] (which has NOT been factorized 
yet!). Thus, one of the major goals in this phase is to predict the required computer memory for 
subsequent numerical factorization. The outputs from this symbolic factorization phase will be 
stored in the following 2 integer arrays (assuming the stiffness matrix data shown in Eq. 2.27 is 


used): 


JSTARTROW 


1 

2 

3 

4 

5 

6 

7 = Y+ 1 


(2.28) 
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JCOLNUM 


(2.29) 


1 \ (4 ’ 

2 6 

3 5 

4 = 1 5 ► 

5 5 

6 6 

7 - NCOEFl j [6 

The following “new” definitions are used in Eqs. (2.28-2.29): 

NCOEF2 s The number of nonzero, off-diagonal terms of the factorized matrix [U] 

JSTARTROW(i) = Starting location of the first nonzero, off-diagonal term for the i“' row of 
the factorized matrix [U]. The dimension for this integer array is N+1 
JCOLNUM(j) s Column numbers associated with each nonzero, off-diagonal terms of [U] 
(in a row-by-row fashion). The dimension for this integer array is 
NCOEF2. Due to “fills-in” effects, NCOEF2 > > NCOEF. 

2.9 SPARSE NUMERICAL FACTORIZATION 

It is generally safe to say that sparse numerical factorization is more complicated for 
computer coding implementation than its skyline, or variable bandwidth cases. Main difficulties 
are due to complex “book-keeping” (or index referring) process. The “key” ideas in this 
numerical phase are still basically quite similar to the operations involved during the symbolic 
factorization phase. There is an important modification that need to be done on the symbolic 
factorization, in order to do the sparse numerical factorization (to facilitate the discussion, please 
refer to the data shown in Eq. 2.27): 

For symbolic factorization purpose, there is no need to have any floating arithmetic 
calculation. Thus, upon completing the symbolic process for row 4, there are practically no 
needs to consider row 2 and/or row 3 for possible contributions to row 5. Only row 4 needs 
to be considered for possible contributions (or “fills-in” effects) to row 5 (since row 4, with 
its “fills-in”, is already full). 

For numerical factorization purpose, however, all rows 2, then 4 and then 3 will have to be 
included in the numerical factorization of row 5. 



36 


2.10 SPARSE SOLVER WITH IMPROVED STRATEGIES 


To simplify the discussion, assuming that upon completion of the symbolic phase, the 
stiffness matrix [K] will have the following form 
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(2.30) 


In Eq. (2.30) the stiffness matrix [K] has 14 dof The symbols “x” and “F” refer to the original 
nonzero terms, and the nonzero terms due to “Fills-in”, respectively. It can be seen that rows 1-3 
have the same nonzero patterns (by referring to the enclosed “rectangular” region, and ignoring the 
fully populated “triangular” region of rows 1-3). Similarly, rows 4-5 have the same nonzero patterns. 
Rows 7-10 have the same nonzero patterns. Finally, rows 11-14 also have the same nonzero 
patterns. Thus, for the data shown in Eq. (2.30), the “Master"’ (or “Super”) dof can be generated as 


MASTER 


1 1 


'3 ' 

2 


0 

3 


0 

4 


2 

5 


0 

6 


1 

7 

zz < 

i 

8 


0 

9 


0 

10 


0 

11 


4 

12 


0 

13 


0 

14 = 


0 


(2.31) 


According to Eq. (2.31), then the “master” (or “super”) dof are dof #1 (which is followed by 2 
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“slave” dof), dof #4 (which is followed by 1 slave dol), dof #6 (which has no slave dof!), dof #7 
(which is followed by 3 slave dof), and dof #1 1 (which is followed by 3 slave dof). 

The skeleton FORTRAN code LDL^ (with full matrix) shown in Table 3 A should be 
modified as shown in Table 3B for incorporating sparse and “Master” (or “Super”) degree-of- 
ff eedom-strategi es 

Table 3B: Pseudo FORTRAN Skeleton Code 
for Sparse LDL^ Factorization With Master dof Strategies 

1. c ... Assuming row 1 has been fa ctorized earlier 

2. DO 1 1 1 = 2, N 

3. DO 22 K = Only those previous “master” rows which have contributions to 
current row 1 

4.1 c... Compute the multiplier(s) (Note: U represent L”^) 

4.2 NSLAVEDOF = MASTER(I) - 1 

5.1 XMULT = U (K,I) / U (K,K) 

5.2 XMUL„, = U (K+m,I) / U(K+m, K+m) 

5.3 c ... m = 1,2, ... NSLAVEDOF 

6. DO 33 J = appropriated column numbers of “master” row #K 

7.1 U (I,J) = U (I,J) - XMULT * U(K,J) 

7.2 - XMUL„, ♦ U (K+m,J) 

8. 33 CONTINUE 

9.1 U (K,I) = XMUL,, 

9.2 U (K+m, I) = XMUL^ 

10. 22 CONTINUE 

11. 11 CONTINUE 

2.11 NUMERICAL EVALUATIONS OF DIFFERENT GENERALIZED EIGEN- 
SOLVERS 

Based upon the discussions in previous sections, practical finite element models (such as 
Exxon off-shore Structure [6], High Speed Civil Transport Aircraft [4], and a 2-D Reverberation 
Chamber [ 1 1 ] are used to evaluate the performance of the developed sparse eigen-solvers. Since 
the codes have been written in standard FORTRAN language (and without using any library 
subroutines), it can be ported to different computer platforms (such as SUN-Sparc-20, IBM- 
R6000/590, Intel Paragon, Cray-C90 etc...) with no (or minimum) changes to the codes. 
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The accuracy of the developed codes for solving generalized eigen equations can be measured 
by the Relative Error Norm (=R.E.N.) which can be computed as: 


R.E.N.= 


1^4 


(32) 


The basic subspace iteration code, given in Ref. [1], will be used as a based-line reference. This 
“basic” subspace iteration code [ 1 ] will be compared to the developed basic, “sparse subspace 
iteration, and “sparse” Lanczos codes. It should be emphasized here that more advanced, faster 
version of Subspace Iteration software has also been reported in the literature [1]. Lumped 
masses are used in all examples in this section. 

Example 1 : Exxon Off-Shore Structure [6] 

The finite element model for the Exxon structure has used extensively in earlier research 
works. The resulted system of generalized eigen-equations from the Exxon model has 23, 155 
dof. 

The number of nonzero terms of the original stiffness matrix is 809,427. Using the Nested- 
Dissection (ND) algorithm, the number of nonzero terms (including “fills’in” terms) is 
10,826,014. The relative error norm (or R.E.N., defined in Eq. 32) and the wall-clock time are 
presented and explained in Figure 2. It should be noted here that on the IBM-R6000/590 
workstation, vector processing capability is available, where as the vector processing capability 
is “not” available on the Sun Sparc-20 workstation. 
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Example 2 : High Speed Civil Transport (HSCT) Aircraft [4] 


The finite element model for the HSCT aircraft has been used extensively in earlier research 
works. The resulted system of generalized eigen equations Irom the HSCT model has 16,1 52 
dof. The number of nonzero terms of the original stiffness matrix is 373,980. Using the 
Modified Minimum Degree (MMD) algorithm, the number of nonzero terms (including “fill-in ’ 
terms) is 2,746,286. The numerical performances of 3 generalized eigen-solvers are presented 
in Figures 3-4. 


Figure 2; Exxon Model 
(USTSU31, Sun-Sparc20 Workstation) 



No. eigen-values 





Figure 3: HSCT Aircraft model 
(Stretch, IBM-R600/590 Worksation) 
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Time (seconds) 


Figure 4: HSCT Aircraft model 
(Rhino, ODU Sun Sparc-20 Worksation ) 
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Example 3 : A Two-Dimensional Reverberation Chamber (R.C.) [11] 


The finite element model for the 2-D Reverberation Chamber (RC) has been used in earlier 
research works [11]. The resulted system of generalized eigen-equations trom the R.C. model 
has 1 0,972 dof The number of nonzero terms of the original “stiffness” matrix is 23 ,946. Using 
the MMD algorithm, the number of nonzero terms (including “fills-in” terms) is 1 65, 135. The 
numerical performances of 3 generalized eigen-solvers are presented in Table 4. 


Table 4: A 2-D R.C. Model 

(Rhino, ODU Sun Sparc-20 Workstation, Consistent mass is used) 


No. Eigen Values 

Norfolk (Sun 3500/ODU), 

Rhino, ODU Sim Sparc-20 

Rhino, ODU Sun Sparc-20 


SV Lan 

SV Sub 

SV Lan 


(in seconds) 

(in seconds) 

(in seconds) 

30 

1 1 .9sec 

221 

51 

100 

170.7 

5,081 

655 

300 

1502.8 (= 25min) 

42,262 

6,696 
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Example 4 : A Simple Mathematical Data For Validating “Complex Numbers” Eigen- 

Solution Case 

In this example one would like to solve the following generalized eigen value ( 1 2 degree-of- 
ffeedom) problem: 

[/f][d)] = U][M][<l)] 


where: 



'( 11 . 0 , 0 . 0 )( 0 ., 0 .) 

( 12 ., 0 .).. 


(o.,o.) 

(13..0.) 

( 22 ., 0 .) 


Jl2xl2 



(l.,0.)(0.,0.) (o.,o.) 

(1..0.) (o.,o.) 

(l.,0.).(0.,0.) 


sym 


( 1 .. 0 .) 


Jl2xl2 


The first 4 “complex” eigen values can be computed from the FORTRAN code as: 

A, = (ll.0,0.0), = Aj = (l3.,0.X andA4 = (l4.,0.) 
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2.12 CONCLUSION 


In this paper, basic generalized eigen-solution algorithms are reviewed. Major 
computational tasks in Subspace iterations, and Lanczos algorithms have been identified. 
Efficient sparse technologies have been developed, and fully utilized, in conjunction with the 
basic Subspace iterations and Lanczos algorithms (such cis: sparse symbolic, numerical 
factorization with unrolling strategies, sparse forward & backward solutions, sparse matrix- 
vector multiplications etc...) for efficient solutions of the generalized eigen-equations. 
Numerical results from practical finite element models have clearly indicated that the 
proposed “sparse” Subspace iterations and Lanczos algorithms have offered substantial 
computational advantages over the traditional “skyline”, or “variable bandwidth” strategies. 

Both lumped and consistent mass formulation have been incorporated into our eigen- 
solution package. Furthermore, both stiffness and mass matrices can be provided with either 
real, or complex numbers. 

Applying these advanced techniques will enable a thorough statistical analysis of the 
reverberation chamber fields. The estimated required order of the eigenvalve problem is very 
large (1 00,000) and is driven by the desire to model the large number of modes which exist at 
higher frequencies. The number of eigensolutions that can be accurately modeled is 
frequency limited. 

Future work will include a complex eigensolver to model the typically lossy behavior of 
reverberation chambers and the resulting statistical field variations. 
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Chapter 3 

Statistical Characterization And The Simulation Of A 
Reverberation Chamber Using Finite Element 
Techniques 

3.1 INTRODUCTION 

A reverberation chamber is an enclosure consisting of metal walls with a metallic paddle 
wheel (denoted a “stirrer” or “tuner”) forming a high quality factor (Q) cavity with continuously 
variable boundary conditions. Reverberation chambers have attained increased importance in the 
determination of electromagnetic susceptibility of avionics equipment. This importance will 
become even more critical as advanced high-speed transport aircraft are developed that 
increasingly depend on electronic sensors and computer control of flight surfaces to manage the 
flight parameters. 

The fields in a reverberation chamber are typically characterized by statistical means. 
Specifically the probability density functions for the real and imaginary components of a 
particular polarization of the electric and magnetic fields are normally distributed. The field 
magnitudes are Rayleigh distributed and the power is exponentially distributed. Reasonable 
statistical agreement for a source-free two-dimensional finite element model has been obtained 
and was the primary emphasis of the work presented by the author [1]. 

There are two focus areas that will be addressed in this work. The first problem to be 
briefly addressed in this paper is an examination of a two-dimensional finite element model for 
both the transverse electric (TE) and transverse magnetic (TM) solutions inside a reverberation 
chamber. A two-dimensional approach to the analysis of reverberation chambers was initially 
suggested by Wu [2] for mechanical stirring using the transmission line matrix (TLM) method. 
Hill [3] examined frequency stirring for an empty two-dimensional structure supporting TM 
modes. For the current work the fields (both TE and TM) in the cavity will be simulated with an 
emphasis on the tuner effects on the modal structure and the resulting statistics of the field 
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distribution. The second focus is shielding effectiveness in a reverberating environment. Typical 
shielding effectiveness measurements are performed in a plane wave environment under various 
angles of incidence. The reverberating environment may provide additional insight into the 
shielding properties in a statistical sense. 

The following section presents an overview of the computational tools used to simulate 
the reverberation environment. The TE results for the simulation “chamber” are then presented 
with an examination of coupling (fields inside vs. outside) shielding effectiveness through an 
aperture. The TM results are then presented for the simulation chamber. The application of a 
two-dimensional analysis tool may lead to useful investigations of stirrer efficiency, field 
homogeneity, and shielding effectiveness. 

3.2 FINITE ELEMENT FUNDAMENTALS FOR TE AND TM FIELDS 

The finite element method is a deterministic approach to the solution of Maxwell’s 
equations using a weighted residual formulation over a set of compact-support basis functions to 
solve for the fields. Consider the geometry of Figure 3.3.1 with a source located at (xo, yo)- 



Figure 3.1. The 2D reverberation geometry 
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The electromagnetic field behavior is governed by Maxwell’s equations as given by 

W y-E = -j<afuH , (31) 

VxH = jcoeE+J. (3-2) 

The inclusion of the source will be accomplished through the characteristics of the electric 
current density, J, in eq. (3.2). By taking the curl of eq. (3.1) and substituting eq. (3.2) and 
assuming non-magnetic media the following inhomogeneous vector wave equation is obtained; 

Vx[VxE)-kle^E = -j\r}o^ (3-3) 

where E is the electric field intensity in volts/meter, is the wavenumber with k^ = , at a 

radian frequency a> , with permittivity s and permeability fi , and intrinsic impedance % . The 
behavior of the electric field is of the form 

E{x,y,z) = \e,{x,y) + ze^ {x, y )] e'''" . 

Setting y = 0 for cutoff and expressing the del operator as V = V, + f , it is possible to wnte 

two separate equations - one for the transverse part and another for the z component. The 
transverse field behavior is modeled with the use of edge elements and the z-directed fields are 
modeled using traditional node based elements. The expression of eq. (3.3) will depend on 
whether a TE or TM field will be considered to exist within the structure. This dependence 

completely rests in the expression of the electric current density, J . 

3.2.1 TE Fundamentals*a dipole source 

The TE case can be expressed by considering a transverse electric current that can be 

thought of as a short dipole antenna with a uniform current density. In the TE case J can be 
expressed 

J = IqJj 

from which the eq. (3.3) is written 

V, X ( V, X ^ = ~jk^T]QlQ Jj (3 -4) 

V,.(V,eJ + *oVz=0- 

Eq. (3.4) represents the transverse variation and the eq. (3.5) represents the z-directed variation. 
The transverse electric current exists only on the edge of a given element and is zero elsewhere. 
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Obtaining the weighted residual form involves integrating over the element support and 
assigning the continuity conditions. Integrating on an elemental basis is accomplished by 

Jj[(V, *1^ ds (3.6) 


and 

Assigning continuity conditions and setting the fields to zero on the perfectly conducting walls 
and tuner yields a matrix equation of the form Ax = b for the transverse fields. Note that for this 

system that x = { } and that Ez = 0. 


3.2.2 TM Fundamentals-an infinite line source 

Following directly the work of Hill [3], the TM case uses an infinite line source located at 

(xo,>'o)ofthe form 


J = zl^S{x-Xo)S{y-yo). 

(3.8) 

so that (4) can be written as 


V,x{V,xe,)-k^£,e, =0 

(3.9) 


(3.10) 

leading to the following weighted residual form: 


Jl[(v, xc,)*(v, =0 

(3.11) 

and 


If [-(V,e, ) • (V,rJ + ys = ) 

(3.12) 


Note that in obtaining eq. (3.12), the integration required to obtain the right-hand side is trivial. 
The resulting linear system of equations is of the form Ax = b . Note that for this system that 

X = and that Ex = 0 and Ey = 0. 


3.3 TE DIPOLE SOURCE RESULTS 

The application of the finite element technique to fields which are transverse electric 
(TE) will be presented in this section. The resulting fields. Ex and Ey , will be used in the two- 
dimensional TE finite element model of the reverberation chamber. Some measures of 
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reverberation chamber performance include the statistical characterization of the fields, stirring 
ratio, field uniformity, the maximum to average ratio, and the normalized standard deviation [4]. 
Several of these measures will be applied to provide basic validation of the reverberating 
characteristics of the simulation “chamber”. The shielding effectiveness in the simulation 
chamber will also be presented in the next sub-section. The 2D dipole source results are obtained 
by the rotation of the tuner for 225 steps in 1.6° increments, thus providing a full mechanical 
rotation of the tuner. The results represent the solution of the matrix equation corresponding to 
eq. (3.6) for each of the tuner positions at each frequency of interest. Consider the geometry of 
Figure 3.3.2 depicting the discretized structure to be analyzed with a typical result in Figure 
3.3.3. Note that the results of Figure 3.3.3 are logarithmic and emphasize the field structure 
rather than the absolute field levels. The source is in the right-hand region at location (2, -1). 



Figure 3.2. Discretized geometry of Figure 3.1. 



'2 X in m, frcq =400 (MHz) 

-2 A 


Figure 3.3. Fields for Figure 3.2 for a fixed tuner 
position. 
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3.3.1 Reverberation in the TE environment 

There are several key factors governing reverberation characteristics and these factors 
can be used to determine if a chamber is operating properly. These factors include the stirring 
ratio, field uniformity, the ratio of the maximum field to average field, and the underlying 
statistical characterization of the fields. Of these measures of reverberation characteristics the 
stirring ratio and field uniformity will not be addressed in this paper. As the values of the field at 
particular point are plotted it was noted that strong numerical resonances occur throughout the 
computational domain in the lossless structure. These resonances correspond to the eigenvalues 
of the source-free problem [1]. In order to improve the numerical stability the quality factor, Q, 
can be altered by filling the chamber with lossy material (e = 1+jO.OOl). The field statistics at a 
given point for 225 tuner steps is depicted in Figure 3.3.4. 



TT3mz rttp pi --0.42, (70,40) 

Wdbull Probability Plot,RdEx2+iniEx2 



Figure 3.4. The statistical characteristics of the total TE field for a 
fixed point for 225 tuner positions. 


The total squared electric field statistics at a point inside the shielded LHS are depicted with a Q 
of 1000 at a frequency of 400 MHz. The Weibull probability plot [1] shown at the bottom of 
Figure 3.3.4 is a strong indicator regarding the nature of a sample statistic. The probability plot 
helps eliminate ambiguities that may arise in interpreting agreement by comparing cumulative 
distribution ftinctions (CDFs). The probability plot maps the parent distribution to a straight line 
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and emphasizes differences in the tails of the sample distribution. Note that the + symbols on the 
plot represent the sample and all but approximately 2 % of the data are nearly collinear with the 
parent Weibull distribution. The Weibull density function is a more general density function that 
can fit a variety of bounded continuous distributions including the exponential distribution. The 
square of a particular component of the electric field, say , can be written 


El =Re{£j' +lm{E^]\ (3-13) 

If the underlying distributions for the real and the imaginary parts are N(0,1) then the squared 
sum is exponentially distributed as well as distributed as chi-squared with two degrees of 
freedom. The probability density function for an exponential distribution is 


T = /W/^) = 



(3.14) 


with parameter //.. The exponential distribution is a special case of a Weibull distribution given 
by the probability density function 




,x>0 
,x < 0 


( 3 . 15 ) 


The exponential distribution corresponds to 6 = 1 in eq. (3.15). Considering the x-component of 
the electric field as depicted in Figure 3.3.4, the estimate of best =1.0799 is a good indicator of 
agreement with the expected distribution. It strongly suggests that the underlying distributions 
are indeed normally distributed. The CDF is plotted in Figure 3.3.5 showing a strong agreement 
with the expected Weibull distribution and the chi-squared distribution with two degrees of 
freedom. Another measure of reverberation quality is the normalized standard deviation cr„ is 
defined by 

<T 

M 

where a is the sample standard deviation and /i is the sample mean. The point shown (70,40) in 
the right hand side (RHS), which is the complex portion of the structure, has a normalized 
standard deviation of 0.991. The ratio of the maximum to average field is 7.93 dB is shown in 
Figure 3.3.6 which agrees well the expected value of 8 dB. The stirrer ratio is 33.8907 dB which 
is indicative of excellent stirring within the complex RHS. A summary of example statistics for 
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a few points in the structure is shown in Table 1. Table 1 portrays the Weibull parameters, the 
normalized standard deviation, the max-to-average ratio, and the stirring ratio for two points in 
each of the two major regions within the simulation chamber. The plots of Figures 4,5, and 6 
represent the data for the point (70,40). Note the general excellent agreement for each of the 
points shown. Of particular interest is that the LHS is unstirred and it has generally identical 
statistical and reverberation characteristics as the complex RHS. This interesting result is 
suggestive of a form of “source” stirring wherein the fields appearing at the aperture joining the 
two regions acts with sufficiently random amplitude and structure as to excite a complex field in 
the non-complex region. 


Table 3.1. Summary of statistical characteristics for the TE fields. 

























Total Field vs. Tuj»r Position , Q *1000, f ■WMHi, Max-to-Avg(dB)*7.9276 



Figure 3.6. The total TE field for a fixed point for 225 tuner 
positions highlighting the max-to-average ratio and stirrer 
ratio. 


Considering the field uniformity within each region can further highlight the idea of 
“source” stirring. In this study the field uniformity will be exploring by considering a test region 
in which the ensemble average of the fields is sampled and compared to other points within the 
test region. The fields will be considered to be uniform if the sample is with ± 3 dB of the 
average of all sample points. The test region is considered uniform if 75% of the sample field 
points are within the + 3 dB region. The total electric field in the RHS is studied in Figure 3.3.7. 
The top of Figure 3.3.7 depicts the sample region wherein 100 points are chosen at random. 
Note that the wall with the aperture is not shown, but exists for reference points (50, Y). The 
uniformity is depicted in the bottom of Figure 3.3.7 with the circled data representing points 
within ± 3 dB of the average value of approximately -7 dB. The triangles indicate points that are 
outside the specified limits. As can be seen, 81 out of 100 points are within the limits and 
therefore the test region in the RHS can be considered to have a uniform field. The uniformity in 
the LHS is shown in Figure 3.3.8 with 91 % of the points falling within ± 3 dB of the average. 
The good result for the field uniformity in the non-complex LHS further suggests that “source” 
stirring can provide excellent reverberation characteristics. 
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Figure 3.7. ±3 (IB Field upiformity pf tbe tQtpl TF field fpr 
100 points ip tite RHS for 225 tpper positions. 




Figure 3.8. ±3 dB Field uniforniity of the total TE field for 
100 points in the XHS for 225 tuner positions. 
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3.3.2 Shielding Effectiveness for the TE structure 

The application to the problem of coupling shielding effectiveness via an aperture will be 
examined in this section. The geometry of Figure 3.1 is the geometry under consideration for 
this examination. Consider a set of random points in the source region (RHS). For those fixed 
points there will be some maximum field level that will exist over all tuner positions. This level 
represents the maximum level in the source region. Consider taking another set of random points 
in the test region (LHS). This set consists of the individual maximums for each test point for all 
tuner positions. A typical result for the shielding effectiveness of the maximums can be seen in 
Figure 3.9. Interestingly there are 22 occurrences of negative shielding effectiveness out of the 
total of 100 random points chosen. This is indicative of phenomena wherein the field is larger in 
the shielded test region than in the source region. The average shielding effectiveness of the 
maximums was 2.19dB for 400 MHz. 


Shielding Effectiveness Max , Q = 1 000, f = 400 MHz 



Figure 3.9. Shielding effectiveness of the maximums for the TE field. 


57 



B 


*■ < 


P 


ly 


3.4 TM INFINITE LINE-SOURCE RESULTS 

The work presented in this section will examine the statistics of the TM fields in the two- 
dimensional reverberation structure. This work closely follows the work of Hill in that the only 
surviving field component is Ez. The results presented in this section will exactly parallel the 
work on the TE fields previously shown, so developmental discussions will be omitted. 


3.4.1 Reverberation in the TM environment 

The field solution for the TM problem consists of the z-component of the electric field 
with real and imaginary parts. The field statistics for the TM case are shown in Figure 3.10 for 
reference point (70,30). The estimate of the shape parameter best = 0.50945 for the Weibull 
distribution indicates that one or more of the underlying distributions is not normally distributed. 
This result is consistent with Hill’s results when the frequency chosen is such that an insufficient 
number of modes exist within the structure. It is significant to note that the field levels for the 
TM problem are several orders of magnitude larger than the TE problem. The CDF of the 
sample is depicted in Figure 3.11 in comparison with the Weibull distribution and the chi- 
squared distribution with the normalized standard deviation of 2.3413. The maximum to average 
ratio is shown in Figure 3.12 and is 12.34 dB. The stirrer ratio is also portrayed in Figure 3.12 
and has a value of 49.56 dB. Table 2 depicts several points in various points within the geometry 
similar to those shown earlier with similar results for the fields in the complex RHS and the non- 
complex LHS. 


RcEz2'^I»Ei3 TV. TomiPofitsm, Q *1000, f*400 MH» 



Figure 3.10. The statistical characteristics of the total 
TM field for a fixed point for 225 tuner positions. 
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R«Es2-»Im£i2 v». Tun*r Pontioa , Q •IfJOO, f «400 IfRx 



Figure 3.11. The cumulative distribution of the total 
TM field for a fixed point for 225 tuner positions. 


Total Field VS. Tuner Position , Q “1000, f Max-to-Avg (dB)"12.3433 



Figure 3.12. The total field total TM field for a fixed point for 
225 tuner positions highlighting the max-to-average ratio and 
stirrer ratio. 
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Table 3.2. Summary of statistical characteristics for the TM fields. 


TM 

Weibull 

Parameters 






Ref. 

Coordinates 

a 

b 

CTn 

MA(dB) 

SR(dB) 


(20,40)- 

(-2.16, -0.42) 

1.042 

0.509 

2.3479 

12.38 

48.68 

} LHS 

(40,40) 

(-0.79, -0.42) 

1.043 

0.508 

2.3461 

12.36 

50.36 

(70,30) 

(1.38,-0.82) 

1.041 

0.509 

2.3413 

12.34 

49.55 

} RHS 

(80,40) 

(2.09, -0.42) 

1.043 

0.508 

2.3454 

12.36 

50.51 


The field uniformity in the RHS and LHS are studied in Figures 13 and 14. Note that the 
RHS has 63 % uniformity in Figure 3.13 and the LHS has 59 % uniformity. These results 
combine to raise some basic questions regarding the statistical characteristics of the TM model at 
the frequencies used, and are the subject of a more detailed exploration in the following 
paragraph. 

Further examination revealed that the real and imaginary parts of is distributed 
according to a Johnson SU distribution [5,6]. The Johnson distribution is denoted Su and is a 
four parameter distribution with a probability density function 
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(3.16) 


for all real numbers x. The parameters are a location parameter y € (— 00 , 00 ) , scale parameter f5 > 
0, and shape parameters a, e (— > and U 2 > 0. The CDF is given by 
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(3.17) 


with 0(z) the standard normal distribution function. The skewness of the density function is 
represented in the shape parameter cc\ with the density function skewed left when cc\ >0, 
symmetrical when a\ = 0, and skewed right when a\ < 0. Tables 3 and 4 depict the summary 
statistics for the real and imaginary parts, respectively. These results depict a wide variation of 
the parameters for various locations within the structure. Work is ongoing to ascertain the 
connection that can possibly be made between these results and an undermoded reverberation 
chamber. 
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Figure 3.13. ±3 dB Field uniformity of the total TM field for 
100 points in the RHS for 225 tuner positions. 




Figure 3.14. ±3 dB Field uniformity of the total TM field for 
100 points in the LHS for 225 tuner positions. 
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Table 3.3. Summary of statistical characteristics for the real part of the TM fields. 


Re {£ 2 } 

Johnson Sy Parameters 


Ref. 

Coordinates 

7 

P 

ai 

«2 

(20,40) 

(-2.16,-0.42) 

14.98 

17.93 

-0.5110 

0.4818 

} LHS 

(40,40) 

(-0.79, -0.42) 

-52.77 

63.63 

0.5249 

0.4903 

(70,30) 

(1.38,-0.82) 

11.93 

26.33 

-0.6526^ 

0.5331 

1 RHS 

(80,40) 

(2.09, -0.42) 

-34.69 

43.59 

0.5366 

0.4931 


Table 3.4. Summary of statistical characteristics for the imaginary part of the TM fields. 


Im{£2} 

Johnson Su Parameters 


Ref. 

Coordinates 

7 

P 

ai 

0.2 

(20,40) 

(-2.16,-0.42) 

-148.64 

199.28 

-0.2355 

0.9527 

} LHS 

(40,40) 

(-0.79, -0.42) 

506.07 

651.37 

0.2439 

0.9361 

(70,30) 

(1.38,-0.82) 

-186.49 

203.49 

-0.2830 

0.8880 1 

} RHS 

(80,40) 

(2.09, -0.42) 

345.93 

^1.17 

0.2482 

0.9342 




3.4.2 Shielding Effectiveness for the TM structure 

The shielding effectiveness of the maximums for the TM field is depicted in Figure 3.15. 
The resulting coupling appears to be much stronger for and identical test configuration. The 
average shielding effectiveness is seen to be 0.991 dB with 40 points out of 100 demonstrating 
negative shielding effectiveness. 


Shielding Effectiveness Max , Q = 1 000, f = 400 MHz 



Figure 3.15. Shielding effectiveness of the maximums for the 
TM field. 
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3.5 CONCLUSIONS 


A two-dimensional finite element model for both the transverse electric (TE) and 
transverse magnetic (TM) solutions inside a simulation reverberation chamber has been 
presented. Tuner effects on the modal structure and the resulting statistics of the field 
distribution have been explored. It was demonstrated that the TE simulation reverberation 
chamber provided excellent statistical and reverberation characteristics. The measures of 
reverberation characteristics included normalized standard deviation, the max-to-average ratio, 
the stirrer ratio, and the field uniformity. In addition the notion of “source” stirring was 
introduced wherein the fields that were coupled into the non-complex geometry were shown to 
be statistically similar to those in the complex environment. The idea of source stirring may 
shed light on the question of whether it is absolutely necessary to stir in all regions that a 
reverberation type test is being performed. The practical issues are enormous when considering 
the demands that internal stirring may place on the testing regimen. The TM simulation chamber 
did not satisfy any of the required reverberation characteristics, and the field components were 
shown to be distributed according to a Johnson Su distribution. Further study may link the 
results of the TM simulation chamber to an undermoded reverberation chamber and may be 
useful in extending the useful low frequency of the reverberation chamber. The shielding 
effectiveness in a reverberating environment was been examined. Coupling shielding 
effectiveness for TE and TM results was presented. The reverberating environment provides 
additional insight into the shielding properties in a statistical sense. 
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Chapter 4 

Shielding Effectiveness In A Reverberation Chamber 
Using Finite Element Techniques 

4.1 INTRODUCTION 

A reverberation chamber is an enclosure consisting of metal walls with a metallic paddle 
wheel (denoted a “stirrer” or “tuner”) forming a high quality factor (Q) cavity with continuously 
variable boundary conditions. Reverberation chambers have attained increased importance in the 
determination of electromagnetic susceptibility of avionics equipment. This importance will 
become even more critical as advanced high-speed transport aircraft are developed that 
increasingly depend on electronic sensors and computer control of flight surfaces to manage the 
flight parameters. 

Shielding of avionics equipment from both interior and exterior electromagnetic threats is 
an important part of the system reliability budget. It is typically assumed that the shielding 
effectiveness of a structure has a level that is a function of frequency for various angles of 
incidence. Measurements in an anechoic chamber and analysis using computational methods 
simulate a plane wave environment to quantify shielding effectiveness. An investigation has 
been performed [1] suggesting that an under certain conditions the aircraft fuselage behaves as a 
reverberation chamber when stirred internally. Another investigation [2] revealed that the fields 
in an aircraft flown in the vicinity of a transmitting antenna have levels with characteristics that 
are statistically similar to a reverberation environment. A conclusion that may be drawn from 
this is that an aircraft fuselage is a structure in which the interior fields are not plane waves. This 
significant connection raises an interesting question: What is the shielding effectiveness of an 
aperture in a reverberating environment? The emphasis of this paper is the examination of a 
two-dimensional finite element model of a reverberation chamber to explore the shielding 
effectiveness of an aperture over a variety of frequencies in a TM field configuration. 
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The fields in a reverberation chamber are typically characterized by statistical means. 
Specifically the probability density functions for the real and imaginary components of a 
particular polarization of the electric and magnetic fields are normally distributed. The field 
magnitudes are Rayleigh distributed and the power is exponentially distributed. Reasonable 
statistical agreement for a two-dimensional finite element model has been obtained and was the 
primary emphasis of the work presented by the author [3] for a source free problem. The field 
statistics for the TM model will not be developed in this work. 

A two-dimensional approach to the analysis of reverberation chambers was initially 
suggested by Wu [4] for mechanical stirring using the transmission line matrix (TLM) method. 
Hill [5] examined frequency stirring for an empty two-dimensional structure supporting TM 
modes by defining the 2D Green’s function for the rectangular structure. For the current work 
the fields in the cavity will be simulated as transverse magnetic (TM) with an emphasis on the 
tuner effects on the modal structure and the resulting shielding effectiveness in a reverberating 
environment. Typical shielding effectiveness measurements are performed in a plane wave 
environment under various angles of incidence. The reverberating environment may provide 
additional insight into the shielding properties in a statistical sense. Park [6] recently published 
results for the simulation of the shielding effectiveness of a two-dimensional TM structure in a 
plane wave environment. Interestingly he compared the simulated two-dimensional results with 
measured results. Based on Park’s work, the shielding effectiveness of the TM structure via 
aperture coupling for several frequencies in a simulated reverberation environment will be 
examined. Measurements were also performed in NASA’s reverberation chamber and are 
presented to provide a connection to the simulated statistical results. 

The following section presents an overview of the computational tool used to simulate 
the reverberation environment. The TM results are then presented for the simulation chamber 
and the replacement (fields in box vs. fields with no box) shielding effectiveness is explored. A 
comparison of the TM results to an anechoic environment is made. Finally a comparison to 
measured results in NASA’s reverberation chamber are presented. These results will point to the 
development of a statistical model of shielding effectiveness that will form a critical component 
to the overall mechanism of upset prediction in digital systems. 
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4.2 FINITE ELEMENT FUNDAMENTALS FOR TM FIELDS 

The finite element method is a deterministic approach to the solution of Maxwell’s 
equations using a weighted residual formulation over a set of compact-support basis functions to 
solve for the fields. Consider the geometry of Figure 4. 1 with a source located at {xo, yo)- 


PECs 


Tuner 



Figure 4.1. Geometry for the TM problem. 

The electromagnetic field behavior is governed by Maxwell’s equations as given by 

VxH = jooeE + J . 

The inclusion of the source will be accomplished through the characteristics of the electric 
current density, J, in eq. (4.2). By taking the curl of eq. (4.2) and substituting eq. (4.1) and 
assuming non-magnetic media the following inhomogeneous vector wave equation is obtained: 

V X ( V X £) - kle,E = 

where £is the electric field intensity in volts/meter, the wavenumber with = ojyJJde , at 
a radian frequency (o , with permittivity s and permeability fu , and intrinsic impedance % . The 
behavior of the electric field is of the form 
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E {x, y, z) = [e, (x, + ze^ (x, )^)] e''" . 

Setting y = 0 for cutoff and expressing the del operator as V = V, +f it is possible to write 

two separate equations - one for the transverse part and another for the z component. The 
transverse field behavior is modeled with the use of edge elements and the z-directed fields are 
modeled using traditional node based elements. The particular form of the expression of eq. 
(4.3) will depend on whether a TE or TM field will be considered to exist within the structure. 
This dependence completely rests in the expression of the electric current density, J . 

4.2.1 TM Fundamentals-an infinite line source 

Following the work of Hill [5], the TM case uses an infinite line source located at 


(xo,yo)of the form 

J = zI^S{x-x^)5{y-y^). (4.4) 

so that (4) can be written as 

V,x(V,xe,)-^„V,=0 (4.5) 

*{'^,e,)^kle^e^=-jk^T]J^5{x-x^S{y-y^) (4.6) 

leading to the following weighted residual form: 

JJ[( V, X e, ) • ( V, X T; ) - ^0 • ^ ] = 0 

and 


Note that in obtaining eq. (4.8), the integration required to obtain the right-hand side is trivial. 
The resulting linear system of equations is of the form Ax = b . Note that for this system that 

X = {E^Y Ex = 0 and Ey = 0. 
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4.3 TM INFINITE LINE-SOURCE RESULTS 

This section presents a TM analysis that will compute the shielding effectiveness of a box 
with a single aperture. The 2D infinite line-source results are obtained by the rotation of the 
tuner for 225 steps in 1 . 6 ° increments, thus providing a full mechanical rotation of the tuner. 
The results represent the solution of the matrix equation corresponding to eq. (4.8) for each of 
the tuner positions at each frequency of interest. Consider the geometry of Figure 4.2 depicting 
the discretized structure to be analyzed with a typical result in Figure 4.3 for a frequency of 400 
MHz. Note that the results of Figure 4.3 are logarithmic and emphasize the field structure rather 
than the absolute field levels. The close-up of the Box in Figure 4.4 demonstrates the field 
structure in the “shielded” box. Note that at 400 MHz the 10 cm aperture is 0.133 wavelengths 
long. Figures 4.5 and 4.6 depict similar results for a frequency of 600 MHz. Note that the field 
structure for 600 MHz is significantly more complex and that the mode structure inside the box 



Figure 4.2. Discretized geometry for Figure 4.1. Grid is 
such that there are 10 samples per wavelength. 

demonstrates the existence of the TMn mode. This indicates that there is one resonance inside 
the box between 400 MHz and 600 MHz, which actually occurs at 480.21 MHz. 
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xoam, Tuner Position -85 


Figure 4.5. TM excitation at 600 MHz of a PEC box with an 
aperture. 




-0 0 -0 7 -0 6 -0 5 -04 -03 -0 2 -0 1 0 0 1 

im m, tuner pW&a -S? 


Figure 4.6. Close-up of PEC box with aperture at 600 
MHz. 
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4.4 SHIELDING EFFECTIVENESS FOR THE TM STRUCTURE 

This section will explore replacement-shielding effectiveness where the fields are 
simulated with and without the box present. This work follows closely the work of Park [6] who 
performed a plane wave analysis of a two-dimensional TM box with single and multiple 
apertures. Park also performed measurements that verified his computational results. This 
section will duplicate Park’s test configuration, but with the significant difference of applying a 
reverberation chamber simulation. Measurements in NASA’s medium size reverberation 
chamber (Chamber B) are presented to further corroborate the simulation. The shielding 
effectiveness (SE) simulation will explore the SE of a two-dimensional box that is 40 cm x 50 

cm with a 10 cm aperture as depicted in Figure 4.1. 

Simulations were performed at various frequencies with the shielding effectiveness 

defined as 


5£ = 20*log 


f NoBox 
r?Box 

\ ^2 J 


(4.9) 


with the NoBox solution corresponding to the solution obtained with no box present and the Box 
solution for the fields at the center of the box. Typical results are portrayed in Figures 4.7 and 


4.8. Figure 4.7 shows a comparison of the z-component of the electric field for the shielded and 
unshielded measurement. Note that these simulated results indicate roughly a 30 dB difference 
between the Box and NoBox cases. The results provide some interesting insight as to the nature 
of SE in a reverberation chamber. Note from Fig. 8 that the SE is normally distributed as seen 
from the linearity of the sample data on the normal probability plot and the subsequent 
histogram. The mean value of the SE is 35.3 dB with a standard deviation of 10.51 dB. Park s 
result is approximately 43 dB and is 7.7 dB greater than the more conservative mean obtained by 
the simulation. The 7.7 dB difference is well within a standard deviation from the simulated 
data. From Fig. 8 it can also be seen that the SE is less than 10 dB for tuner positions 25 and 71 , 
and is nearly zero for position 7 1 . The variation in SE is not attributable to resonance conditions 
since the first mode of the box is at 480.23 MHz. The fields in the aperture are varying 
continuously in amplitude and structure as the tuner is stirring the fields. Although the simulated 
frequency is below the cutoff of single mode operation in the box a statistical distribution of the 
field is obtained inside due to what may be considered “source” stirring. This form of stirring is 
distinct from frequency stirring [4] since the frequency in the simulation is fixed. The fields in 
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the aperture have a statistical distribution governed by the characteristics of the simulation 
reverberation chamber and provide an excitation to the smaller non-complex cavity. This 
excitation drives the non-complex box in a manner consistent with large, complex reverberation 
chambers. An important conclusion to the idea of source stirring is that internal stirring of the 
non-complex box is not required. 


U 


|E^| vs. Tuner Position,Frequency =400 MHz 



Figure 4.7. Magnitude of for the shielded and unshielded cases at 
400 MHz at the center of the box. 
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4.4.1 Measurements 

This section will highlight the measurement of shielding effectiveness in the 
reverberation chamber. These measurements were performed on a rectangular box (40 cm x 50 
cm X 20 cm) with an aperture that is 10 cm x 2 cm as depicted in Figure 4.9. 



Figure 4.9. Box with aperture for shielding effectiveness 
measurements in the reverberation chamber. 


For this measurement an HP 8720C vector network analyzer (VNA) was used to measure the 
transmission characteristics (52i) in the chamber. The VNA was calibrated using “Response and 
Isolation”. The measurement consists of a log-periodic antenna (transmit) to a probe (receive) 
located in the box at three positions: top, side, and rear wall for shielded measurements and then 
a measurement without the box for an unshielded measurement. The SE can be computed by a 
slight variation of eq. (4.9). 


20* log 


f o NoBox \ 
^2\ 

^Box 

V ^^21 J 


(4.10) 


for each wall position. A typical result at 400 MHz is depicted in Figure 4.10. A comparison of 
the SE reported by Park [6] and the simulated and measured results are presented in Figure 4.11. 
The Park data overestimates the SE for all but 500 MHz and is 28 dB higher than the simulated 
finite element model at 100 MHz. The simulated and measured results are nearly identical at 
300 and 400 MHz, but diverge at the low and higher frequencies. The simulated results are 17 


76 


dB higher than the measured results at 100 MHz, and are 5 dB lower at 600 MHz. Further 
analysis would readily discern the resonance observed in the measured data at 450 MHz. 



- 100 ' ' ' ' ' ' 

0 50 100 150 200 

Tuner step 

Figure 4.10. IS 21 I at 400 MHz for the shielded and unshielded configurations. 
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Figure 4.11. A comparison of reported [6J, simulation chamber, and measured 
shielding effectiveness. 
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4.5 CONCLUSIONS 


A two-dimensional finite element model for transverse magnetic (TM) solution inside a 
reverberation chamber has been presented. Tuner effects on the modal structure and the 
shielding effectiveness in a reverberating environment has been examined. The reverberating 
environment provides additional insight into the shielding properties in a statistical sense. The 
shielding effectiveness of the TM structure for aperture coupling for several frequencies in a 
simulated reverberation environment was presented and compared to published and measured 
results. These results point to the development of a statistical model of shielding effectiveness. 
This model will form a critical component to the overall mechanism of upset prediction in digital 
systems. 

An example of an application of the notion of a probabilistic model for shielding 
effectiveness a problem of determining the probability of exceeding a particular voltage for a 
given input voltage. Consider the problem of determining the probability of having an electric 
field intensity of 1 V/m anywhere inside a cavity exposed at 1000 V/m. The shielding 
effectiveness was characterized earlier as a Normal random variable with mean of 35.3 dB and 
standard deviation of 10.5 dB at a frequency of 400 MHz. The probability of exceeding IV/m 
for a 1000 V/m exposure corresponds to calculating the probability of having a shielding 

effectiveness less than 60 dB which is Fj£^(60) = l — 0l — — 1 = 0.9906. This calculation 

predicts a 99.06% probability of exceeding IV/m for an exposure of 1000 V/m. If the threat 
level is reduced to 100 V/m the probability of exceeding 1 V/m reduces to 67.36 %. 
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